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Ph. We investigate numerical solution of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi 

(DGLAP) Q 2 evolution equation for the transversity distribution A T q or the structure 
function hi. The leading-order (LO) and next-to-leading-order (NLO) evolution equa- 
tions are studied. The renormalization scheme is MS or MS in the NLO case. Dividing 
the variables x and Q 2 into small steps, we solve the integrodifferential equation by 
the Euler method in the variable Q 2 and by the Simpson method in the variable x. 
Numerical results indicate that accuracy is better than 1% in the region < x < 0.8 
if more than fifty Q 2 steps and more than five hundred x steps are taken. We provide 
a FORTRAN program for the Q 2 evolution and devolution of the transversity distri- 
bution A T q or hi. Using the program, we show the LO and NLO evolution results of 
the valence-quark distribution A T u v + A T d v , the singlet distribution ^2 i (A T q i + A T cji), 
and the flavor asymmetric distribution A T u — A T d. They are also compared with the 
longitudinal evolution results. 
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Program Summary 



Title of program: H1EVOL 

Computer: AlphaServer 2100 4/200 (SUN-IPX); Installation: The Research Center for 
Nuclear Physics in Osaka (Saga University Computer Center) 

Operating system: Open VMS V6.1 (SunOS 4.1.3) 

Programming language used: FORTRAN 77 

Peripherals used: Laser printer 

No. of lines in distributed program: 1203 

Keywords: Polarized parton distribution, transversity distribution, chiral-odd structure 
function, Q 2 evolution, numerical solution. 

Nature of physical problem 

This program solves the DGLAP Q 2 evolution equation with or without next-to- 
leading-order (NLO) a s effects for a transversely polarized parton distribution, so called 
transversity distribution A T q or chiral-odd structure function h±. 

Method of solution 

The DGLAP integrodifferential equation is solved by dividing the variables x and Q 2 
into very small steps. The integrodifferential equation is solved step by step by the 
Euler method ( "brute- force" method) in the variable Q 2 and by the Simpson method 
in the variable x. 

Restrictions of the program 

This program is used for calculating the Q 2 evolution of various transversity distri- 
butions A T q (hi). The Q 2 evolution equation is the DGLAP equation. The double 
precision arithmetic is used. The minimal subtraction (MS) or the modified minimal 
subtraction (MS) scheme is used in the NLO evolution. The NLO evolution is identical 
in both schemes. A user provides the initial distribution as a subroutine or as a data 
file. Examples are explained in Section f|. Then, the user inputs seventeen parameters 
in Section [|. 

Typical running time 

Approximately forty seconds on AlphaServer 2100 4/200 and two and a half minutes 
on SUN-IPX. 
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LONG WRITE-UP 



1 Introduction 

The European Muon Collaboration (EMC) shed light on the proton spin structure by 
measuring the longitudinally polarized structure function g-y. Its conclusion seemed 
startling: almost none of the proton spin is carried by quarks. Since then, many 
investigations have been done for understanding the internal spin structure of the 
proton. However, we should wait for refined experiments in order to determine each 
parton polarization, in particular polarized sea-quark and gluon distributions |J. 

In addition to the abovementioned studies on the longitudinal polarization, it is 
interesting to test the spin structure in the transverse polarization. There exists a 
leading-twist structure function which will be measured in the transversely polarized 
Drell-Yan processes |2|, |[ |J. It is now named as the h\ structure function ||, which 
measures the transversity distribution, A T q or often denoted as Sq. It has a chiral-odd 
property so that it cannot be found in inclusive deep inelastic electron scattering. The 
hi structure function is expected to be measured in the Relativistic Heavy Ion Collider 
(RHIC) Spin project || at Brookhaven National Laboratory. The experimental results 
will provide us a good opportunity to understand unexplored transverse spin physics at 
high energies. Furthermore, semi-inclusive electron deep inelastic scattering could also 
reveal the hi although there is a complication of unknown fragmentation functions |J . 
With these experimental possibilities in mind, theoretical physicists should investigate 
detailed properties of hi |7[] as far as they can before the completion of the RHIC 
facility. 

There are, in general, two variables in the structure functions: Q 2 and the momen- 
tum fraction x. In the transversity distribution measured by the Drell-Yan, Q 2 is the 
dimuon-mass square: Q 2 = m 2 . The x distribution is associated with nonperturbative 
physics, which can be handled only in certain quark models at this stage because there 
is no experimental information on h\. The lattice QCD calculation is not still very 
reliable as far as the x distribution is concerned ||. On the other hand, the Q 2 de- 
pendence can be predicted in perturbative Quantum Chromodynamics (QCD) unless 
Q 2 is small. This scaling violation is described by the Dokshitzer-Gribov-Lipatov- 
Altarelli-Parisi (DGLAP) equation ||, which is an integrodifferential equation. The 
leading-order (LO) DGLAP equation for hi was derived in Ref. [O and its numerical 



analysis is discussed in Ref. In addition, the next-to-leading order (NLO) analy- 
sis for its anomalous dimensions and splitting function was completed recently for the 
transversity distribution A T q. The NLO evolution of A T q was studied in the Feynman 
gauge (T^j and in the light-cone gauge [jl3| . It was also investigated in the Feynman 



gauge in a subsequent paper |TJ], |TJ]. The results made it possible to investigate the 
details of NLO evolution effects in the structure function hi [[U| [17], fL8|. In particular, 
the NLO analysis is usually used in the unpolarized structure functions and in the 
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longitudinally polarized structure function g\. The NLO studies can be extended now 
to the hi case. In this paper, we report our studies on the numerical solution. 

Our group has been investigating the numerical solution in the unpolarized structure 
functions and also in the polarized structure function g\. So far, we employed the 
Laguerre-polynomial and the brute-force p(J methods. Other methods with the 
Mellin transformation are, for example, discussed in Ref. f2lfl . The variables Q 2 and 
x are divided into small steps, then the integrodifferential equation is solved step by 
step in Ref. ||20|| . Because the scaling violation is a small logarithmic effect, the 
integration over Q 2 is accurately done with a small number of steps. In this paper, we 
change the brute-force integration over x to the Simpson method with expectation of 
better accuracy. The transversity evolution is rather simple in the sense that the gluon 
does not couple to the distribution A T q because of the chiral-odd nature. Therefore, 
the evolution equation is a single integrodifferential equation without mixing with the 
gluon distribution. It is the same form as the nonsinglet evolution equation. This fact 
simplifies our numerical analysis. The transversity evolution is completed by solving 
only the "nonsinglet-type" evolution equation. 

The purpose of our study is to provide a useful Q 2 evolution program for the 
distribution A T q (or hi). It should be very useful for theoretical and experimental 
researchers who are involved in the RHIC-Spin type project or in semi-inclusive lepton 
scattering one for finding hi. We explain the details of our studies in the following. 
In Section |2j, the transversity DGLAP Q 2 evolution equation is explained. Then, 
our numerical-solution method is described in Section |^. Input parameters and input 
distributions are discussed in Section |4], and subroutines in the program are explained 
in Section [5]. Numerical results and their comparison with longitudinally polarized 
distributions are discussed in Section || Summary is given in Section |7[ Explicit forms 
of splitting functions are listed in Appendices. 
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2 Q 2 evolution equation for the transversity 
distribution A T q 

The transversity distribution A T qi, which is equivalently denoted as h\ or Sqi, with the 
quark flavor % can be measured, for example, in the double transverse-spin asymmetry 
Att in the Drell-Yan process 0, [| |J: 

sin 2 # cos2(j)^2 i e 2 i A T q i (x 1 ) A T qi(x 2 ) . . 

^TT — — j— : oTj 2 ? A - ? \ • \ ZA > 

l + cos 2 6» }^e; f ft(xi) gi(x 2 ) 

The factor (Xrprp — sin # cos 20/(1 + cos 2 0) is the parton's double asymmetry The 
angles 9 and <fi are the polar and azimuthal angles of the lepton momentum with 
respect to the beam and proton polarization directions respectively According to the 
above equation, the distribution combination A T qi A T (ji weighted by the charge square 
can be measured with information on the unpolarized distributions qi and q\. The spin 
asymmetry Att will be measured in the RHIC-Spin project ||, and the transversity 
distributions could be obtained by Eq. fl2.1|) . 

The distribution A T q is interpreted in a parton model: it is the probability to find 
a quark with spin polarized along the transverse spin of a polarized proton minus the 
probability to find it polarized oppositely (A T g = —qi). The LO evolution equation 
for A T q was derived in Ref. [fT0[| , and the NLO evolution of A T q was recently obtained in 



Refs. [[12], [13], |14 |. Because of the chiral-odd nature of the transversity distribution A T q 
or the structure function hi, the gluon does not participate in the evolution equation. 
Therefore, the DGLAP equation is a single integrodifferential equation: 



9 a ±/ _ a s{Q 



2\ 



. nii02 A T g ± (x,Q i ) = ^-±A T P q± {x)® A T g ± (x,Q 2 ) , (2.2) 

where A T P q ± is the splitting function for the transversity distribution. The convolution 
<g> in Eq. fl2T2]) is defined by 

^ dz „ fx 



f( X )®g( X )= / _fl_\g( Z ) . (2.3) 



The notation g ± in the splitting function indicates the A T q + = A T q + A T q or A T q~ = 
A T q — A T q distribution type. The a s (Q 2 ) is the running coupling constant. Both 
LO and NLO evolution calculations can be handled by Eq. ( |2.2|) ; however, there are 
differences between the LO and NLO splitting functions and also between the coupling 



constants. The LO and NLO expressions for A T P g ± [0 and a s are listed in Appendix 



A. In the following numerical analysis, we use the variable t which is defined by 

t = lnQ 2 (2.4) 
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instead of Q 2 . Furthermore, the parton distribution and the splitting function multi- 
plied by x, 

f(x) = xf(x) , (2.5) 
are used throughout this paper because they satisfy the same integrodifferential equa- 



tion. Then, Eq. (2.2) becomes 



A T q ± (x,t) = ^A T P q± (x)®A T q ± (x,t) . (2.6) 



3 Simpson method 



There are various numerical methods for solving the DGLAP equations [21]. We have 



investigated two methods, the Laguerre-polynomial [plj and the brute-force |2(J meth- 



ods. The splitting functions and parton distributions are expanded by the Laguerre 
polynomials in the first method. Even though the Laguerre method has an advantage 
of computing time, it may not be very accurate in the spin-independent nonsinglet 
case at small x. In the brute-force method, the variables t and x are divided into small 
steps, then the integral is calculated in the simplest manner. This method is accu- 
rate if a large number of steps is taken. However, it inevitably takes much computing 



time. In order to resolve these problems, the Simpson method |22[ is employed in the 
z integration of Eq. ( |2.3| ). Here, we report its numerical analysis. 

In the following numerical solution, the brute-force method is still used in the 
variable t and the Simpson method is introduced in calculating the z integration. The 
brute-force method for solving a differential equation is often called Euler method [p2]| , 
so that we use this terminology in this paper. The variables t and z are divided into 
N t and N z steps, then differentiation and integration are defined by 



df(t) _ f(t j+1 ) - /(*,-) 



dt 6t 



(3.1) 



g(z) dz=^2 — [g(z k -i) + 4#(z fc ) + g(z k+1 

k=2 6 



where k = 2, 4, • • •, N z . (3.2) 

The k summation is taken over the even numbers from 2 to N z , so that N z has to be 
an even number. 

With these replacements, the evolution equation can be solved rather easily. Equa- 
tion (|2.6|) is written in the following form: 

A T q (x m ,t j+1 ) =A T q (x m ,tj) + dt— > —F(z k ,x m ,tj) 

271 * — ' 6Zk 

k=2m 

k = 2m, 2m + 2, • • •, N z , (3.3) 
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where the variable x is divided into N x = N z /2 steps. The function F is defined by 
the following: 



F(z k ,x m ,tj) = A T P q ±(z k _ 1 )A T q 



( 



Zk-l 



m 



, tj ) + 4 A T P q ±(z k ) A T q ± pi ^ 



+ A T P g ±(2; fe+1 ) A T g 




(3.4) 



First, the evolution from t\ = lnQp to t% — t\ + St is calculated in the above equation 
by providing the initial distribution A T q ± (x m / Zk,tj =1 ). Repeating this step N t — 1 
times, we obtain the final distribution at t = tjv t +i- However, the small x region is 
often important practically. If the accurate results are required in the small x region, 
x and z should be replaced by log 10 x and log 10 z. Then, the evolution equation Eq. 
(|2.6|) is replaced by 



In order to apply the Simpson method to this equation, log 10 x and log 10 z are divided 
into N x and N z steps. Our FORTRAN code supports both the linear-step and the 
logarithmic-step cases, and the details are explained in the following sections. 



0_ 

Ot 



A T q ± (log 10 x,t) = lnlO 



a s (t) 



f 

Jlo 



d(1og w z)A T P q ±(\og w z)A T q ± (log 10 -,t) .(3.5) 
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4 Description of input parameters and initial 
distribution 



For running the FORTRAN-77 program H1EV0L, a user should supply seventeen 
input parameters from the file #10. In addition, an input distribution(s) should be 
given in a function subroutine(s) in the end of the FORTRAN program or in an input 
data file(s). The i-th initial distribution is written in the output file #30+2 (#31 - 
#38). Evolution results for the distribution-? are written in the output file #20+i 
(#21 - #28). The maximum number of the input distributions is eight. We explain 
the input parameters and the initial distributions in the following. 



4.1 Input parameters 

There are following seventeen parameters. 



IREP 


1 


finish the program after this evolution 




2 


proceed to the next evolution with new input parameters 


IOUT 


1 


write x and xA T q(x) at fixed Q 2 (=Q2) in the file(s) #21 - #28 




2 


write Q' 2 and xA T q(Q 2 ) at fixed x (=XX) in the file(s) #21 - #28 


IREAD 


1 


give initial distribution(s) in function subroutine(s) 




2 


give initial distribution(s) from the file(s) #11 - #18 


INDIST 


1 


do not write initial distribution(s) 




2 


write initial distribution(s) in the file(s) #31 - #38 




3 


write initial distribution(s) in the file(s) #31 - #38 
without calculating evolution 


IORDER 


1 


leading order (LO) in a s 




2 


next-to-leading order (NLO) 


IMORP 


1 


xA T qi — xA T q~i type distribution 




2 


xA T qi + xA T qi type distribution 


ILOG 


1 


linear- a; and linear- .z steps are taken in calculating the evolution 
and in writing the x-dependent output. 




2 


logarithmic- a; and logarithmic- z steps 

(We recommend to use ILOG=2 in a general case.) 



initial Q 2 (= Qq in GeV 2 ) at which an initial distribution is supplied 



Q02 



Q' 2 to which the distribution is evolved (Q 2 # Qq, Q' 2 > Q'q or Q 2 < Q'fy 
QCD scale parameter Aqcd in GeV 



Q2 



DLAM 



number of quark flavors Nf 



NF 



x at which Q' 2 dependent distributions are written (IOUT=2 case) 



XX 



number of x steps N x (NX < 3000) 



NX 



NT 



number of t steps N t (NT < 3000) 



NSTEP 



number of x steps or t steps for writing output distribution(s) 



XMIN 



minimum of x, x min (0 < XMIN < XX) 



NFI 



number of distributions which are evolved simultaneously (NFI < 8) 
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Numerical values of the parameters should be supplied in the file #10, then these are 
read in the program through the subroutine GETPA1. It should be noted that the first 
four lines of the file #10 are comments. Therefore, the parameters are given from the 
fifth line as shown in the following example. There are irrelevant parameters in some 
cases: 

IOUT = 1 case ■ ■ ■ XX 
INDIST = 3 case ••• I ORDER, IMORP, Q2, DLAM, NF, XX, NT 
IORDER = 1 case ■ ■ ■ IMORP . 

Although they are not used in running the evolution part, their numerical values should 
be supplied in the file #10 within the allowed ranges. 

The meaning of IREAD is explained in Section 4.2. The parameter IMORP indi- 
cates a plus or minus type distribution ^2 i aix(A T q i ± A T qi), where Oj are some con- 
stants. In the NFI > 2 case, distributions with the same ± type can be evolved at the 
same time. The maximum number of the distributions is set up as eight in the program. 
For example, if one would like to evolve two A T q + type initial distributions at Q 2 =4.0 
GeV 2 to the distributions at Q 2 =200.0 GeV 2 by the NLO DGLAP equation with N f =A 
and A=0.231 GeV, the input parameters could be IREP=1, IOUT=l, IREAD=1, 
INDIST=1, IORDER=2, IMORP=2, ILOG=2, Q02=4.0, Q2=200.0, DLAM=0.231, 
NF=4, XX=0.1, NX=500, NT=50, NSTEP=50, XMIN=0.0001, and NFI=2. In this 
case, the input file #10 is: 



IREP,IOUT,IREAD,INDIST,IORDER,IMORP,ILOG 
Q02,Q2,DLAM,NF,XX,NX,NT,NSTEP,XMIN,NFI 



1, 1, 1, 1,2, 2,2 

4.0, 200.0, 0.231, 4, 0.1, 500, 50, 50, 0.0001, 2 

The first four lines are comments. The parameter values from IREP to ILOG should 
be written in the fifth line and the remaining ones in the next line. If one would like 
to repeat the evolution, for example if the LO evolution is also needed, the input file 
is the following: 



IREP,IOUT,IREAD,INDIST,IORDER,IMORP,ILOG 
Q02,Q2,DLAM,NF,XX,NX,NT,NSTEP,XMIN,NFI 



2, 1, 1, 1,2, 2,2 

4.0, 200.0, 0.231, 4, 0.1, 500, 50, 50, 0.0001, 2 
1, 1, 1, 1, 1,2,2 

4.0, 200.0, 0.231, 4, 0.1, 500, 50, 50, 0.0001, 2 
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In the same way, one can repeat the evolution further by choosing IREP=2. The 
evolved results are written continuously in the same output file(s). 



4.2 Initial distributions supplied by function subroutines 
(IREAD=1) 

If IREAD=1 is chosen, an input distribution(s) at Ql should be supplied in the end of 
the FORTRAN program H1EVOL as a function subroutine(s) HlINl(X) - H1IN8(X). 
For example, three distributions (e.g. xA T u + , xA T d + , and xA T s + ) can be evolved at 
the same time in the NFI=3 case. These initial distributions should be supplied in 
the function subroutines HlINl(X), H1IN2(X), and H1IN3(X). The remaining func- 
tions H1IN4(X) - H1IN8(X) may be set to zero. The evolution results for the initial 
distribution-i, which is given by the function HlINi(X) (i=l, 2, • • •, NFI), are written 
in the output file #20+i. If INDIST=2 or 3 is chosen, the initial distribution(s) is 
written in the file #30+i. 

To be exact, there is no information on the input transversity distribution because 
experimental data do not exist. We discuss this point in Section 6.1. In our program, 
the Gehrmann-Stirling (GS) type-A [Q singlet distribution xA T q s = (z^-r^i + 
xA T q,j) is given in the function HlINl(X), and the distribution [(xA T u + xA T u) — 
(xA T d + xA T d)}/2 is given in H1IN2(X). Additionally, the valence-quark distribution 
xA T u v + xA T d v is given in the function HlINl(X) and the distribution [(xA T u — 
xA T u) — (xA T d — xA T d)}/2 is given in H1IN2(X) as comments. 



4.3 Initial distributions supplied by data files 
(IREAD=2) 

If IREAD=2 is chosen, an input distribution(s) at Ql should be supplied in a separate 
data file(s) # 11 - #18 as shown in the following. We give the singlet distribution as 
an example. 



0.000100 
0.000110 
0.000120 
0.000132 
0.000145 



0.001721 
0.001612 
0.001485 
0.001339 
0.001172 



1.000000 



0.000000 



The first column is the x values and the second one is the corresponding distribution 
values. The data should be sorted from the smallest x to larger ones, and the data 
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at x < Xmi n and at x=1.0 must be supplied. The last line number should be smaller 
than three thousand. The evolution results for the initial distribution in the file #10+i 
(i=l, 2, • • NFI) are written in the file #20+i. If INDIST=2 or 3 is chosen, the initial 
distribution(s) at NSTEP+1 points is written in the file(s) #30+i. 

5 Description of the program H1EVOL 

5.1 Main program H1EVOL 

The seventeen input parameters are read by calling the subroutine GETPAl and these 
parameters are checked by the subroutine ERR1. If there is an error, the program 
writes the error message in the output file #6 and stops. Color constants and other 
necessary constants are calculated in the subroutine GETPA2. Then, t, x, z, and the 
splitting function are evaluated at each step point by the subroutines GETT, GETX, 
GETZ, and GETP respectively. If IOUT=2 is chosen, the XMIN is replaced by XX 
and x values in the array X(I) are shifted in the subroutine GETMIN so as to have 
X(1)=XX, X(NX+1) = 1.0. The initial distribution(s) is calculated in the subroutine 
GETINI or is read from the data file(s) by the subroutine REDINI, then it is stored in 
the array HI (IN, I). Here, the argument IN is the distribution number from 1 to NFI. 
The initial distribution(s) is written in the file(s) #30+IN by the subroutine OUTXD 
if INDISiyi. 

Since all the necessary values are ready in the above procedures, we start the evo- 
lution calculation. First, the initial distribution(s) is stored in the array WH1(IN,I) 
and a s at t\ is calculated by the subroutine GETALP. The evolved distribution(s) at 
t 2 — t\ + St is calculated by the subroutines GETINT and GETH1, and the results 
are stored in HI (IN, I). If IOUT=2, the result (s) at x — XX is stored in the array 
H1Q2(IN,I). Next, the evolution results are stored in WH1(IN,I) and a s at t 2 is eval- 
uated. Then, the evolution from t 2 to £3 = t 2 + St is calculated in the same way. 
Repeating this step, we obtain the final distribution(s) at t — ijv t +i- The calculated 
values are interpolated either in t or in x depending on IOUT by using the subroutine 
SPLINE and the function SEVAL. The final results at NSTEP+1 points are written 
in the output file(s) #10+IN by the subroutines OUTXD or OUTQ2D. 

5.2 Subroutines GETPAl (IREP,IEORD,IERR), GETPA2, 
and ERR1 (IREP,IERR) 

The subroutine GETPAl reads the seventeen input parameters from the input file #10. 
These parameters are checked whether they are within the valid ranges by the subrou- 
tine ERR1. If there is an error, IERR=1 is returned from the subroutine ERR1 and 
the program stops with warning notification in the file #6. The subroutine GETPA2 
calculates necessary color constants, n, {3q, and (5\. 
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5.3 Subroutines GETT(Q02,Q2,NT,T,DT) 

and GETX(XMIN,XMAX,NX,ILOG,X,DX) 

These subroutines evaluate the values of t = In Q 2 and x at N t + 1 and N x + 1 points, 
and they are stored in the arrays T(I) with 1=1, 2, • • •, N t + 1 and X(I) with 1=1, 2, • • •, 
N x + 1. In the subroutine GETX, either the linear-a; (IL0G=1 case) or logarithmic- 
x (ILOG=2 case) step is chosen. These subroutines also return the value of St = 
[t{Q 2 ) - t(Ql)]/N t = DT and 5x = DX, where 5x is defined by Sx = (1 - x min )/N x in 
the ILOG=l case or 5x = (0 — log 10 x m j„)/A^ in the ILOG=2 case. 

5.4 Subroutine GETMIN(XX,X,NX,XMIN) 

In the IOUT=2 case, the subroutine GETMIN is called after getting X(I). This subrou- 
tine replaces the XMIN value by XX because the distribution information at x < XX 
is not necessary in the evolution calculation. Then, the values in the array X(I) are 
shifted. For example, if XX is in the region X(M) < XX < X(M+1), X(l) is replaced 
by XX and X(2), X(3), • • •, X(NX-M+2) are replaced by X(M+1), X(M+2), • • •, 
X(NX+1). The number of steps N x is also replaced by N x - M + 1. 

5.5 Subroutine GETZ(ILOG,X,NX,DX,Z,NZ,DZ) 

The subroutine GETZ calculates the values of z at N z + 1 (N z = 2N X ) points and 
stores these values into the array Z(K), K=l, 2, • • •, N z + 1. The first one is given by 
Z(l) = X(l), and the second one Z(2) is calculated by Z(2) = Z(l) + [X(2) -X(l)]/2 in 
ILOG=l case or Z(2) = Z(l) + [log 10 X(2) - log 10 X(l)]/2 in ILOG=2 case. The sub- 
sequent ones are given by Z(3)=X(2), Z(A)=Z(3)+5x/2, Z(5)=X(3), Z(6)=Z(5)+5a;/2, 
• • • in the ILOG=l case or by Z(3)=X(2), log 10 Z(4)=log 10 Z(3)+5a;/2, Z(5)=X(3), 
log 10 Z(6)=log 10 Z(5)+<5x/2, • • • in the ILOG=2 case. 

5.6 Subroutine GETP(IORDER,SIGN) 

The subroutine GETP calculates the LO and NLO splitting functions at the N z + 1 
points and stores the calculated values in the arrays P0(K), POPF(K), P1(K) and 
P1PF(K). The splitting functions are calculated by the functions PLO, POPL, PNLO, 
and P1PL. Additionally, the extra factors associated with the '+' function or 5 function 
are also evaluated by calling the subroutines PODL and P1DL. 

5.7 Subroutines GETINI(NFI,H1) and REDINI(NFI,H1,IERR) 

If IREAD=1 is chosen, the subroutine GETINI calculates the initial distribution(s) 
by calling the function(s) HlINl(X) - H1IN8(X). The subroutine REDINI reads the 
initial distribution(s) from the input data file(s) and interpolates the distributions(s) 
by using SPLINE and SEVAL if IREAD=2. 
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5.8 Subroutine GETALP(T,IORDER,ALP2PI) 

The subroutine GETALP calculates the a s /27r=ALP2PI at given t. 

5.9 Subroutines GETINT(DINTEG), FUN(K,I,IN), 
and SIMPS(FUN,I,IN,DZ) 

The subroutine GETINT calculates the integral by calling the subroutine SIMPS. The 
SIMPS integrates the integrand FUN by using the Simpson method. The integral 
results in the array DINTEG are then returned. The factor IN is the input distribution 
number, IN=1, 2, • • -, or NFL The factors K and I represent the z and x positions, 
Z(K) and X(I). 

5.10 Subroutine GETH1(DT1, DINTEG, Hll) 

The subroutine GETH1 calculates the distribution at t = tj from the one at t — tj-i 
with the integral obtained by GETINT. The evolution results are stored in Hll (IN, I). 
The factor DTI is the t step. 

5.11 Subroutines OUTXD(ILOG,Hl,NFILE) 
and OUTQ2D(Q2T,HlQ2) 

These subroutines interpolate the final evolution results and write the distribution(s) 
at NSTEP+1 points in the output file(s). The OUTXD (OUTQ2D) writes the x (Q 2 ) 
dependent results. The OUTXD is also used for writing the initial distribution(s) in the 
INDIST=2 or 3 case. The NFILE is the factor for writing either the initial distribution 
(NFILE=1) or the final one (NFILE=2). 

5.12 Subroutine SPLINE(N,X,Y,B,C,D,IN) 
and function SEVAL(N,XX,X,Y,B,C,D,IN) 

The subroutine SPLINE calculates the coefficients, B(IN,I), C(IN,I), and D(IN,I) (I = 
1, 2, • • •, N) by the cubic Spline interpolation for the distribution type IN. The inter- 
polated points Y(IN,I) at X(I) (/ = 1, 2, • • •, N) should be supplied. This interpolation 
program is taken from Ref. p4| . Using the obtained Spline coefficients, the function 
SEVAL calculates the value of Y at given XX. 

5.13 Functions PLO(Z), POPL(Z), PNLO(Z,SIGN), 
and P1PL(Z) 

The functions PLO and PNLO calculate the LO and NLO splitting functions. The 
functions POPL and P1PL calculate the LO and NLO '+' function parts in the splitting 



12 



functions. If IMORP=l (SIGN=-1.0), the PNLO calculates the A T q - A T q type 
splitting function, and it calculates the A T q + A T q type if IMORP=2 (SIGN=+1.0). 



5.14 Subroutines P0DL(PLUS0,DEL0) and P1DL(PLUS1,DEL1) 

The subroutines PODL and P1DL calculate the extra factors which are related to the 
'+' function [see the second term of Eq. ( A.8|) ] and the 5 function. These factors are 
taken out from the integral. 



5.15 Functions S2(Z), SPENCE(X), and DGAUSS(EXTERN,A,B) 

The function S2 calculates S^x) which appears in the splitting function. The ^(x) 
is expressed in terms of the Spence function in Eq. ( |A.12j ). The Spence function is 
defined by the integral form in Eq. ( |A.13|) , and it is numerically calculated by the 
function DGAUSS in the Gauss-Legendre method with the integrand SPENCE. 



5.16 Function GS(X,A,B,C,D,E,F) 

This function calculates the GS type parton distribution with the parameters A, B, C, 
D, E, and F as ABx c {l - x) D (l + Ex + Fy/x). 



5.17 Functions HlINl(X), H1IN2(X), • • •, and H1IN8(X) 

These functions calculate the initial distributions. For example, three initial distribu- 
tions should be provided in the NFI=3 case by HlINl(X), H1IN2(X), and H1IN3(X), 
and the remaining functions should be set to zero. 
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6 Numerical analysis 



6.1 Accuracy of Q 2 evolution results 

There are two parameters which determine the accuracy of our numerical results. They 
are the numbers of steps: N x and N t . Another number of steps N z is fixed by the 
relation N z = 2N X . In this subsection, we discuss how the evolution results depend on 
these parameters and how long it takes for calculating the evolution. 

There is a practical problem in the transversity evolution in the sense that there 
is no input distribution. For example, the input parton distributions are obtained by 
fitting various experimental data in the unpolarized distribution studies. At this stage, 
there is no experimental information on the transversity distribution, so that the initial 
distribution could not be supplied. An item of good news is that the lattice QCD studies 
started producing a x distribution || although it is still somewhat away from the very 
reliable one. Therefore, we may use the initial distribution calculated by a bag model 
a nonrelativistic quark model |16|, 2l|, or the Gliick-Reya-Stratmann- Vogelsang 



distribution [26, 11 at very small Q 2 . It is known that the transversity distribution 
and the corresponding longitudinally polarized one are almost the same at small Q 2 
according to Ref. ||. However, the applicability of perturbative QCD becomes dubious 
in the small Q 2 region. Therefore, instead of stepping into the nonperturbative region, 
we decided to use a longitudinally polarized distribution at larger Q 2 as the input 
transversity distribution. The Gehrmann-Stirling (GS) type-A distributions [23] are 
used throughout this paper. They are given at Q 2 =4 GeV 2 . We warn the reader that 
the assumption Aq = A T q would not be valid in the Q 2 =4 GeV 2 region even though it 
would work at small Q 2 (~ Aq CD ). Future experimental data will clarify the difference 
between the two distributions. 

We show the dependence on N x and N t by taking the singlet distribution. First, we 
fix the x-step number at A^=1000, then the t-step number N t is varied as 10, 50, 200, 
and 1000. The input distribution is the GS-A singlet distribution at Q 2 =4 GeV 2 [p3}[ . 
Because there is little information on each antiquark distribution ||27|| , flavor symmetric 
distributions are assumed for the antiquark distributions. Calculated evolution results 
for the distribution xA T q s = x^2 i f (A T qi + A T qi) are shown in Fig. [I]. We use 
Nf =4 and A=231 MeV in calculating the evolution. The GS-A singlet distribution is 
evolved to the one at Q 2 =200 GeV 2 by our NLO program. The initial distribution is 
supplied at the end of the program as the subroutine HlINl(X). The input parameters 
are IREP=1, IOUT=l, IREAD=1, INDIST=1, IORDER=2, IMORP=2, ILOG=2, 
Q02=4.0, Q2=200.0, DLAM=0.231, NF=4, XX=0.1, NX=1000, NT=10, 50, 200, or 
1000, NSTEP=100, XMIN=0.0001, and NFI=1. There is almost no difference between 
the various results in Fig. |l] except for the A^=10 case. It means that merely fifty steps 
are enough for the getting the accurate evolution. The small number of steps is enough 
because the scaling violation is a small logarithmic effect. 

Next, the step N x is varied with fixed ^=100 in Fig. 0, where ^=50, 200, 500, and 
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3000 are taken. Numerical accuracy is slightly worse in the N x =50 case; however, the 
accurate results are obtained by taking the five hundred steps. Defining the evolution 
accuracy by \[A T q(N x ,N t ) - A T q(N x = 3000, N t = 1000)]/ A T q(N x = 3000, N t = 
1000) |, we find that the accuracy is better than 1% in the region 10~ 5 < x < 0.8 with 
N x = 500 and Nt = 50 in the logarithmic- a; steps (ILOG=2). If one is more interested 
in the large-re region, one had better take the linear- a; steps (ILOG=l). However, we 
recommend the user to use the logarithmic steps in a general case. 

Finally, we discuss a typical running CPU time. If one initial distribution is evolved 
with N x = 500 and Nt = 50, it is about forty seconds on the AlphaServer 2100 
4/200 and two and a half minutes on the SUN-IPX. From these analyses, we find that 
our program can be run on a typical workstation or possibly on a personal computer 
without spending much computing time. However, if one would like to calculate the 
evolution repeatedly or to evolve several distributions simultaneously, a reasonably 
powerful machine is recommended. 



6.2 Q 2 evolution of transversity distributions 



There are some studies on the transversity evolution fllT , IB|, [IT], IB| ; however, it is not 
a well investigated topic. Therefore, we show evolution results of various transversity 
distributions in this subsection. In particular, the LO and NLO evolution differences are 
shown, then they are compared with those of the corresponding longitudinally polarized 
distributions by using the BFP1 program ||20|| . The transversity NLO evolution is the 



same in the MS |T2| , 14 and MS |H| schemes. 



First, the singlet evolution results are shown in Fig. ^. The same GS-A distribu- 
tion is assumed for the transversity and longitudinally-polarized parton distributions 
at Q 2 =4 GeV 2 . Furthermore, the LO and NLO distributions are assumed the same 
in both cases. The initial distribution is shown by the dotted curve. It is evolved to 
the distributions at Q 2 =200 GeV 2 by the transverse or longitudinal evolution equa- 
tion. The transversity NLO effects increase the evolved distribution at medium-large 
x and also at small x (<0.01), and they decrease the distribution in the intermediate x 
region (0.01 < x < 0.1). The NLO contributions are rather different from those of the 
longitudinal evolution. The transversity evolution differs significantly from the longi- 
tudinal one either in the LO case or in the NLO. The evolved transversity distribution 
A T q s is significantly smaller than the longitudinal one Aq s in the region x ~ 0.1. The 
magnitude of A T q s itself is also smaller than that of Aq s at very small x (< 0.07). 
Therefore, one should be careful about the Q 2 evolution difference between the two 
distributions in extracting information on A T q s from experimental data. The explicit 
Q 2 dependence is shown in Fig. f|, where the distributions are shown at £=0.001, 0.01, 
0.1, and 0.5. 

In parametrizing the parton distributions [|IJ, it is natural to provide the valence and 
sea quark distributions separately. We expect that future parametrizations on A T q s 
will be given in the same manner. Therefore, it is important to show the evolution 
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results of a valence-quark distribution. The evolution calculations are done in the 
same way as those in Fig. [3] except that the initial distribution is the nonsinglet one 
x(Au v +Ad v ) — x(A T u v +A T d v ). The evolution results are shown in Fig. [|. The dotted 
curve is the initial distribution, and evolved distributions are solid and dashed curves. 
The evolved transversity distributions are significantly smaller than the longitudinally 
polarized ones, in particular at small x. 

The other interesting distribution is the flavor asymmetric distribution A T u — A T d. 



Although the unpolarized version u — d is rather well known [|27] due to the mea- 
surements of the Gottfried-sum-rule violation and the Drell-Yan p-n asymmetry, the 
longitudinally polarized one Au — Ad is not investigated experimentally. Because of 
the difference between the splitting functions for the A T q + A T q and A T q — A T q type 
distributions, a finite A T u — A T d could be obtained through the Q 2 evolution [111] as it 



appears also in the unpolarized case ]27J. We show the evolution results of A T u — A T d 
and Am — Ad in Fig. ^| As it is obvious from the figure, the obtained distributions 
are fairly small. In order to get the accurate evolution in Fig. |6], large numbers of 
x and t steps should be chosen: for example iV x .=3000 and iV 4 =500. There is no fla- 
vor asymmetry in the GS-A antiquark distributions, so that the initial distributions 
are A T u — A T d = Au — Ad = 0. Because the difference between the q ± splitting 
functions does not appear in the LO, the evolved distributions show no asymmetry 
(A T u — A t (1)lo = (Au — A(1)lo = at Q 2 =200 GeV 2 . The situation is completely 
different in the NLO. The NLO evolution results are shown by the solid and dashed 
curves. It is interesting to find the finite distributions in the NLO evolution. The 
evolved transversity distribution x(A T u — A T d) is concentrated in the x ~ 0.04 region, 
and the longitudinal one is distributed in the smaller x (~0.004) region. However, 
we should note the magnitude. The distributions are typically 0.0001, which is too 
small to be found experimentally in the near future. The situation is the same in the 
unpolarized: the perturbative QCD correction to the Gottfried sum is typically 0.3% 
with respect to the 30% experimental violation. The A T u — A T d as well as Am — Ad 
will be measured in the W-charge-asymmetry experiments at RHIC [P^| . Even if fi- 
nite distributions are found, they are unlikely to be explained by the perturbative 
mechanism. 
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7 Summary 



We investigated numerical solution of the Q 2 evolution equation for the transversity 
distribution A T q or the structure function hi. A useful FORTRAN program was 
created for calculating the evolution in the LO and NLO cases. The renormalization 
scheme is MS or MS in the NLO. The numerical solution was obtained by dividing 
the variables x and t = In Q 2 into small steps with the numbers of steps, N x and N t . 
The Euler and Simpson methods are used for solving the integrodifferential equation. 
We find that the evolution accuracy is better than 1% with N x > 500 and N t > 50 in 
the x region 10~ 5 < x < 0.8 . 

Using the program, we showed the evolution results of the valence-quark distri- 
bution x(A T u v + A T d v ), the singlet distribution xj^i / (^ T 5'« + A T gj), an d the flavor 
asymmetric distribution x(A T u — A T d). The evolution results are compared with those 
calculated by the longitudinal DGLAP evolution equations. Because the transverse and 
longitudinal evolution results are significantly different, we should be careful in ana- 
lyzing experimental data. The Q 2 variations of the flavor asymmetric distribution are 
fairly small so that a nonperturbative mechanism should be explored if a finite distribu- 
tion is found experimentally. The transversity evolution of these distributions should 
be tested by future experimental projects such as the RHIC-Spin. 
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Appendix A. Splitting functions and running 
coupling constants 

The running coupling constant in the leading order (LO) is given by 



47T 



/3 hi(Q 2 /A 2 ) 

and the one in the next-to-leading order (NLO) is 



(A.l) 



a 



NLO 



(Q 2 



47T 



A)hi(Q 2 /A 2 



l 



Aln{ln(Q 2 /A 2 )} 
/? 2 ln(Q 2 /A 2 ) 



(A.2) 



The constants (3q and (5\ are defined by the color constants Cq, Cp, and Tr as 

A = y C g - \t r N s , A = ^C 2 G - ™C G N f - 2C F N f , (A.3) 

where 



a 



G 



c h 



2N C 



R 



(A.4) 



with the number of color (N c =3) and the number of flavor (Nf). 

The splitting function is expressed by the LO and NLO splitting functions as 



A T P q± (x) = A T P^(x) + ^^A T P ( 



(A.5) 



The second term is the NLO contribution. We use the expressions in Ref. |L3) for the 
transversity splitting functions. The LO splitting function is given by |L0| 



A T P£\x) = C F 

where the '+' function is defined by 

/(*) 



2x 



x 



+ ^S(l-x) 



(A.6) 



dz ■ 



:i-z) 



1 — z 



(A.7) 



It should be noted that the above integration is defined in the region < z < 1. 
However, the integration is given in the region x < z < 1 in the actual evolution 
calculation. In such a case, the integral with '+' function is calculated by the following 
equation: 



(1 - z)- 



l dz f ^ + /(l)ln(l-x) 

1 — z 



(A.8) 
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The NLO splitting function is [13 



A T P£>(x) = A T pU(x) ± A T P 



<i<i 



qq \ 



(A.9) 



1 — x 



+ 21n(l-:z)J \nxS T P^>{x) 



+ 



3 - 2 



7T 



+ 6C(3) 5(1 -x) 



67 11 



+ ( ^ + ^ In x + In 2 x - ^ ) 5 T P( Q) ( x ) + ( + 



9 



+ -C F T R N f 



:i-x) 



fn n 



7T 



V 12 



9 



6C(3)U(l-x) 



lnx-^5 T Pi g 0) (x)-Q + ^-)5(l-x) 



A r P«(a;) = CV(CV-^C< G 



:1-.t)+25 2 (x)5 t p4 0) (^./-) 



(A.10) 
(A.ll) 



where 5xPqq\x) = 2x/(l — x) + . The C, function is defined by ((k) = Y^n=x^-/ nk -> 
and C(3) is given by the numerical value (£(3)=1. 2020569...). The ± factor in Eq. 
(|A.9| ) indicates "A T g ± A T q type" distribution Yli a i^T < ii ^ ^-t?*)- The ^2 function 
is expressed in terms of the Spence function S'(x) as 



_ , . f 1 ^ dz 1 - 
^(x) = / — In 

Jx/(l+a>) z : 

= 5 f ^— ] - S 
\l + x 

where S(x) is defined by 



x 



In 2 



1 + x 



In 2 



x 



1 + x 



(A.12) 



1 , \nz 

dz 

1 - z 



(A.13) 



The NLO splitting function or anomalous dimensions are the same in the MS scheme 



1^, [n] and in the MS scheme |13|]. Therefore, our evolution program can be used in 
both schemes. 



19 



References 



[1] RHIC-SPIN-J working group on parametrization: Y. Goto, N. Hayashi, M. Hirai, H. 
Horikawa S. Kumano, M. Miyama, T. Morii, N. Saito, T.-A. Shibata, E. Taniguchi, 
and T. Yamanishi, research in progress. 

[2] J. P. Ralston and D. E. Soper, Nucl. Phys. B152 (1979) 109. 

[3] R. L. Jaffe and X. Ji, Phys. Rev. Lett. 67 (1991) 552; Nucl. Phys. B375 (1992) 527; 
X. Ji, Phys. Lett. B284 (1992) 137. 

[4] J. L. Cortes, B. Pire, and J. P. Ralston, Z. Phys. C55 (1992) 409; J. C. Collins, S. 
F. Heppelmann, and G. A. Ladinsky, Nucl. Phys. B420 (1994) 565; B. Kamal, Phys. 
Rev. D53 (1996) 1142; R. L. Jaffe and N. Saito, Phys. Lett. B382 (1996) 165. 

[5] Proposal on Spin Physics Using the RHIC Polarized Collider (RHIC-Spin collabo- 
ration), August 1992; update, Sept. 2, 1993. 

[6] P. J. Mulders and R. D. Tangerman, Nucl. Phys. B461 (1996) 197. 

[7] B. L. Ioffe and A. Khodjamirian, Phys. Rev. D51 (1995) 3373; J. Soffer, Phys. 
Rev. Lett. 74 (1995) 1292; S. Aoki, M. Doui, T. Hatsuda, and Y. Kuramashi, Phys. 
Rev. D56 (1997) 433; A. Y. Umnikov, H. He, and F. C. Khanna, Phys. Lett. B398 
(1997) 6; R. Kirschner, L. Mankiewicz, A. Schafer, and L. Szymanowski, Z. Phys. 
C74 (1997) 501; M. Meyer-Hermann and A. Schafer, |hep-ph/9709349t K. Suzuki and 



T. Shigetani, [hep-ph/9709"39l R. L. Jaffe, X. Jin, and J. Tang, |hep-ph/9709"322 



[8] M. Gockeler, R. Horsley, H. Perlt, P. Rakow, G. Schierholz, A. Schiller, and P. 
Stephenson, [hep-ph/9711245| . 

[9] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15 (1972) 438 and 675; G. 
Altarelli and G. Parisi, Nucl. Phys. B 126 (1977) 298; Yu. L. Dokshitzer, Sov. Phys. 
JETP 46 (1977) 641. 

[10] X. Artru and M. Mekhfi, Z. Phys. C45 (1990) 669. 

[11] V. Barone, T. Calarco, and A. Drago, Phys. Lett. B390 (1997) 287. 

[12] S. Kumano and M. Miyama, Phys. Rev. D56 (1997) 2504. 

[13] W. Vogelsang, |hep-ph/9706511| , Phys. Rev. D in press. 

[14] A. Hayashigaki, Y. Kanazawa, and Y. Koike, |hep-ph / 9707208| , Phys. Rev. D in 
press; numerical analysis in progress. 

[15] The NLO evolution was first reported in Ref. [12| as the Los Alamos preprint |hep- 



ph/9706420| . The first version had some typing and calculation mistakes. The correct 
final result was reported in the subsequent preprint |hep-ph / 9707208! Jl4| within the 
same formalism as well as in Ref. |12| . 

[16] S. Scopetta and V. Vento, |hep-ph/9707250 . 



20 



[17] C. Bourrely, J. Soffer, and 0. V. Teryaev, |hep-ph/9710224 . 



[18] O. Martin, A. Schafer, M. Stratmann, and W. Vogelsang, [hep-ph/9710300 



[19] G. P. Ramsey, J. Comput. Phys. 60 (1985) 97; J. Bliimlein, G. Ingelman, M. Klein, 
and R. Riickl, Z. Phys. C 45 (1990) 501; S. Kumano and J. T. Londergan, Comput. 
Phys. Commun. 69 (1992) 373; R. Kobayashi, M. Konuma, and S. Kumano, Comput. 
Phys. Commun. 86 (1995) 264. 

[20] M. Miyama and S. Kumano, Comput. Phys. Commun. 94 (1996) 185; M. Hirai, 
S. Kumano, and M. Miyama, |hep-ph /9 70 7220| , Comput. Phys. Commun. in press. 

[21] J. Bliimlein, B. Geyer, and D. Robaschik, |hep-ph/9711i05| . 

[22] B. Carnahan, H. A. Luther, and J. O. Wilkes, Applied Numerical Methods (John 
Wiley & Sons, 1969). 

[23] T. Gehrmann and W. J. Stirling, Phys. Rev. D. 53 (1996) 6100. 

[24] G. E. Forsythe, M. A. Malcolm, and C. B. Moler, Computer Methods for Mathe- 
matical Computations (Prentice-Hall, 1977). 

[25] B.-Q. Ma, I. Schmidt, and J. Soffer, |hep-ph/ 9710247 . 

[26] M. Gliick, E. Reya, M. Stratmann, and W. Vogelsang, Phys. Rev. D53 (1996) 
4775. 

[27] S. Kumano, |hep-ph/ 9702367 . 

[28] N. Saito, pp. 40-49 in Spin Structure of the Nucleon, edited by T.-A. Shibata, 
S. Ohta, and N. Saito, (World Scientific, 1996); J. Soffer and J.-M. Virey, |hep 
ph/9706229 |; B. Kamal, |hep-ph/ 97 10374 . 



21 



TEST RUN OUTPUT 



IOUT= 1 IRE AD = 1 INDIST= 1 IORDER= 2 IMORP= 2 ILOG= 2 

Q02= 4.0000 Q2= 200.000 DLAM= 0.2310 NF= 4 NFI= 1 

XX= 0.1000000 NX= 500 NT= 50 NSTEP= 50 XMIN= 0.0001000 
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Figures 




Figure 1: N t dependence of evolution results is shown for the singlet transversity 
distribution. The initial distribution is assumed to be the longitudinally polarized 
GS-A distribution at Q 2 =4 GeV 2 . It is evolved to the one at Q 2 =200 GeV 2 by the 
next-to-leading-order DGLAP evolution equation with Aqc£>=231 MeV and Nf—A. 
Aa.=1000 is fixed and N t is varied as 10, 50, 200, and 1000. There are dotted, dashed, 
dot-dashed, and solid curves at Q 2 =200 GeV 2 for A t =10, 50, 200, and 1000 respectively. 
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Figure 2: N x dependence of evolution results is shown for the singlet transversity 
distribution. The initial distribution is the same as the one in Fig. [I]. iV t =i00 is fixed 
and N x is varied as 50, 200, 500, and 1000. There are dotted, dashed, dot-dashed, and 
solid curves at Q 2 =200 GeV 2 for A^=50, 200, 500, and 1000 respectively. 



25 



0.25 




0.0001 0.001 0.01 0.1 1 

X 



Figure 3: The singlet transversity evolution results are compared with those of longi- 
tudinally polarized singlet distribution. The longitudinal evolution is calculated by the 
BFP1 program pOfl . The same GS-A distribution is assumed at Q 2 =4 GeV 2 , and it is 
shown by the dotted curve. It is evolved to the distributions at Q 2 =200 GeV 2 . The 
dashed and solid curves are the LO and NLO evolution results respectively. 
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Figure 4: Q 2 dependent results of the transversity distribution are compared with those 
of the longitudinally polarized one. They are shown at x=0.001, 0.01, 0.1, and 0.5. The 
solid (dashed) curves are the NLO (LO) evolution results of xA T q s . The dot-dashed 
(dotted) curves are the NLO (LO) evolution results of xAq s . 
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Figure 5: Evolution results of the valence-quark transversity distribution x(A T u v + 
A T d v ) are compared with those of the longitudinally polarized one. The initial GS- 
A distribution at Q 2 =4 GeV 2 is shown by the dotted curve, and it is evolved to the 
distributions at Q 2 =200 GeV 2 . The dashed and solid curves are the LO and NLO 
evolution results respectively. 
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Figure 6: Evolution results of the flavor asymmetric transversity distribution x(A T u — 
A T d) are compared with those of the longitudinally polarized one. There is no flavor 
asymmetry in the initial GS-A distribution at Q 2 =4 GeV 2 [x(A T u — A T d)= x(Au — 
Ad)=0]. In the LO evolution, the evolved distributions are identically zero [x(A T u — 
A T d) LO = x(Au - Ad) LO =0]. The NLO evolution results at Q 2 =200 GeV 2 are shown 
by the solid and dashed curves for the x(A T u — A T d) and x(Au — Ad) distributions 
respectively. 
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